{
  "metadata": {
    "kernelspec": {
      "name": "xpython",
      "display_name": "Python 3.13 (XPython)",
      "language": "python"
    },
    "language_info": {
      "name": ""
    }
  },
  "nbformat_minor": 5,
  "nbformat": 4,
  "cells": [
    {
      "id": "3fbc5e0c-5d55-440a-8de9-e4dda817361a",
      "cell_type": "code",
      "source": "import numpy as np\nfrom qutip import Qobj, sigmax, sigmay, sigmaz, mesolve, Options\nimport scipy.optimize\n\n# Define physical constants\nh = 6.626e-34  # Planck's constant (J·s)\nc = 3.00e8     # Speed of light (m/s)\nhbar = h / (2 * np.pi)\n\n# Parameter ranges\nomegas = np.array([0.05e6, 0.1e6, 0.2e6, 0.5e6])  # Rabi frequencies (MHz)\ngammas = np.array([0.05e6, 0.1e6, 0.2e6])          # Decoherence rates (MHz)\nwavelengths = np.array([630e-9])                   # Single wavelength for simplicity (m)\nomega_0 = 2 * np.pi * c / wavelengths[0]           # Transition frequency (Hz)\n\n# TLS Hamiltonian\nsigma_x, sigma_y, sigma_z = sigmax(), sigmay(), sigmaz()\nH0 = 0.5 * omega_0 * sigma_z\nH1 = 0.5 * sigmax()\n\n# Time parameters\ntlist = np.linspace(0, 50e-6, 500)  # Time array (µs)\n\n# Function to solve dynamics and extract steady-state\ndef get_steady_state(omega, gamma):\n    H = H0 + 0.5 * omega * H1\n    psi0 = Qobj([[1], [0]])\n    c_ops = [np.sqrt(g) * s for g, s in zip([gamma], [sigma_x, sigma_y, sigma_z])]\n    result = mesolve(H, psi0, tlist, c_ops=c_ops, e_ops=[sigmaz()], options=Options(store_states=True))\n    steady_pop = 0.5 * (1 - result.expect[0][-1])  # Last time point\n    return steady_pop\n\n# Function to fit coherence decay (exponential decay model)\ndef exp_decay(t, A, tau):\n    return A * np.exp(-t / tau)\n\ndef fit_coherence_decay(omega, gamma):\n    H = H0 + 0.5 * omega * H1\n    psi0 = Qobj([[1], [0]])\n    c_ops = [np.sqrt(g) * s for g, s in zip([gamma], [sigma_x, sigma_y, sigma_z])]\n    result = mesolve(H, psi0, tlist, c_ops=c_ops, e_ops=[sigmax(), sigmay()], options=Options(store_states=True))\n    coherence = np.sqrt(result.expect[0]**2 + result.expect[1]**2)\n    popt, _ = scipy.optimize.curve_fit(exp_decay, tlist, coherence, p0=[1.0, 1e-6])\n    return popt[0], popt[1]  # Amplitude A, decay time tau (s)\n\n# Generate derived datasets\nsteady_state_data = []\ncoherence_fit_data = []\nsample_excerpt_data = []\n\nfor omega in omegas:\n    for gamma in gammas:\n        # Steady-state population\n        steady_pop = get_steady_state(omega, gamma)\n        steady_state_data.append([omega * 1e-6, gamma * 1e-6, steady_pop])\n        \n        # Coherence decay fit parameters\n        A, tau = fit_coherence_decay(omega, gamma)\n        coherence_fit_data.append([omega * 1e-6, gamma * 1e-6, A, tau * 1e6])  # tau in µs\n        \n        # Sample excerpt (first 10 time points)\n        H = H0 + 0.5 * omega * H1\n        psi0 = Qobj([[1], [0]])\n        c_ops = [np.sqrt(g) * s for g, s in zip([gamma], [sigma_x, sigma_y, sigma_z])]\n        result = mesolve(H, psi0, tlist[:10], c_ops=c_ops, e_ops=[sigmaz()])\n        excited_pop = 0.5 * (1 - result.expect[0])\n        for t, pop in zip(tlist[:10] * 1e6, excited_pop):  # Time in µs\n            sample_excerpt_data.append([t, pop, omega * 1e-6, gamma * 1e-6])\n\n# Save to CSV (example for steady-state, adapt for others)\nimport pandas as pd\nsteady_df = pd.DataFrame(steady_state_data, columns=['Omega (MHz)', 'Gamma (MHz)', 'Steady_State_Excited_Population'])\nsteady_df.to_csv('tls_steady_state_populations.csv', index=False)",
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "execution_count": null
    },
    {
      "id": "bf5f325d-d5c3-48d6-8b29-07a93a0310ff",
      "cell_type": "code",
      "source": "",
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "execution_count": null
    }
  ]
}